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Abstract 

We present a complete next-to-leading order (NLO) calculation for the total cross section of 
inclusive Higgs pair production via bottom-quark fusion (bb — ► hh) at the CERN Large Hadron 
Collider (LHC) in the Standard Model. The NLO QCD corrections lead to less dependence on 
the renormalization scale (/j,r) and the factorization scale (hf) than the leading-order (LO) cross 
section, and they significantly increase the LO cross section. The rate for inclusive Higgs pair 
production is small in the Standard Model, but can be large in models with enhanced couplings of 
the b quark to the Higgs bosons. 
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I. INTRODUCTION 



In the Standard Model (SM), one Higgs doublet is responsible for the electroweak sym- 
metry breaking (EWSB) that generates masses for gauge bosons and fermions. A neutral 
Higgs boson (h) remains after EWSB, and it is the only SM elementary particle that has 
not been observed in high energy experiments. In extensions of the Standard Model, there 
can be more Higgs bosons. 

One of the most important goals of the Fermilab Tevatron Run II and the CERN Large 
Hadron Collider (LHC) is to discover the Higgs bosons or to prove their non-existence. The 
present lower bound on the standard Higgs boson mass from direct searches at LEP2 0, 0] 
is Mh > 114 GeV. The electroweak precision measurements set an upper limit of Mh < 
166 GeV at 95% confidence level for the Standard Model Higgs boson [3| using the recently 
measured top quark mass of m t = 171.4 ± 1.2 ± 1.8 GeV This limit increases to 199 GeV 
when the LEP2 direct search limit is included. 

The Fermilab Tevatron and the LHC will play crucial roles in Higgs searches. Once a 
candidate Higgs boson is discovered, it will be necessary to determine the Higgs couplings 
and spin to see if the Higgs candidate has the properties of the SM Higgs boson. One of the 
most difficult properties to measure is the trilinear self coupling of a Higgs boson [E IH IE HI ■ 
The high energy and high luminosity at the LHC might provide opportunities to detect a 
pair of Higgs bosons as well as a discovery channel to measure the trilinear Higgs couplings 
in the SM and in models with more Higgs bosons. In the Standard Model, gluon fusion 
is the dominant process to produce a pair of Higgs bosons via triangle and box diagrams 
with internal top quarks and bottom quarks h], 11 , 0, ^| . Bottom quark fusion can also 
produce Higgs pairs at a lower rate. 

At tree level, the physical production mechanism for a Higgs boson pair in association 
with b quarks is gg — > bbhh. This process contains a large collinear logarithm from the gluon 
splitting into a collinear bb pair, A = \n(Mh/mb). These logarithms can be resummed by 
using a perturbatively defined b quark parton distribution function (PDF) which is inherently 
0(a s A) [IE IE IEE3- I n this approach, the ordering of perturbation theory is changed to 
be an expansion in 0{pt s ) and A -1 . 

When using a scheme with b quark PDFs for Higgs pair production in association with 
b quarks, the leading order (LO) process becomes bb — > hh, and we compute the NLO cross 
section with 0(a s ) and 0(1/A) corrections to this process. The subprocess bg — » bhh is a 
correction of 0(1/ A) to the lowest order process, while gg — > bbhh is 0(1/A 2 ). 

The rate for Higgs pair production at the LHC is small in the SM. However, it can become 
significant in models in which the Higgs coupling to the bottom quark is enhanced [18]. In 
two Higgs doublet models with Model II type of Yukawa interactions (including the minimal 
super symmetric standard model (MSSM)), the ratio of the Higgs vacuum expectation values 
(tan f3 = v 2/Vl) is an important parameter. A large value of tan f3 greatly enhances the Higgs 
coupling with bottom quarks and makes bottom quark fusion become the largest production 
mechanism for producing a Higgs boson at the LHC. The DO experiment has placed limits 
on single Higgs production in association with b quarks for large values of tan/3 figj ] . 

In this paper, we present the complete next-to-leading order (NLO) calculations to the 
production of Higgs pairs via bottom quark fusion in the Standard Model. We compute a 
consistent set of 0(a s ) and 0(1/ A) corrections. In a future paper, we will present results 
for double Higgs production at NLO in the MSSM 



The theoretical prediction depends on the number of b quarks tagged. Here, we consider 
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only inclusive processes in which there are no tagged b quarks. We apply the two cutoff 
phase space slicing method 21, 2^| to calculate corrections from real gluon emission. Two 
arbitrary small parameters are introduced to split the phase space in real gluon emission 
into soft, hard/collinear, and hard/non-collinear regions. The production in the hard/non- 
collinear region is finite and can be calculated numerically. Divergences are isolated into 
soft and hard/collinear regions. The soft (infrared) singularities cancel with the infrared 
singularities in the virtual corrections. The collinear singularities are absorbed into the 
initial parton distribution functions. 

In section II, we compute the leading-order cross section for pp — > hh + X via bb — > 
hh. In section III, we provide a complete next-leading order (NLO) calculation for hh 
production. At the parton level, we compute one-loop virtual corrections and real gluon 
emission corrections to bb — > hh and the 0(1/ A) subleading process, bg — > bhh. The 
associated production, gg — » bbhh, is finite and is a subleading correction of 0(1/ A 2 ) to the 
inclusive rate for pp — > hh + X. We use MadGraph |23| and HELAS 24] to compute its 
LO production rate with a finite mass of the bottom quark. Numerical results are given 
in Section IV and conclusions are drawn in section V. In addition, there is an appendix to 
present formulas for the b quark running mass. 



II. LEADING-ORDER CROSS SECTION FOR bb -> hh 



The leading order (LO) inclusive cross section for pp — > hh + X via bb —+ hh is evaluated 
with 

cr LO (s,t,u)(bb -> hh) (1) 

where b(x) and b(x) are the LO parton distribution functions for bottom quarks in the proton, 
&Lo( s ,t, u ) is the parton level cross section for bb — > hh and s,t,u are the Mandelstam 
variables. Fig. Q shows the tree level Feynman diagrams for bb — > hh. 




FIG. 1: The lowest order Feynman diagrams for bb — > hh. 



dxidx 2 b(xi)b(x 2 ) + b(xi)b(x 2 ) 



We assign momenta to the initial and the final state particles with 

Kpi)KP2) h(P3)h(p4.) 

and pi + p 2 = Ps + Pa- We take the bbh and hhh couplings to be —V^-C^h and —iS-^Chhh, 
respectively, with v = the Higgs vacuum expectation value ~ 246 GeV, Cbbh = 1 and 
Chhh = 1 in the Standard Model. We evaluate the bottom quark mass in the hbb Yukawa 
coupling by using the MS mass, fn b (fj,) (defined in Appendix A), for a two-loop heavy 
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quark running mass [25| with m^pole) = 4.7 GeV and the NLO evolution of the strong 
coupling [26l | . 

In addition,the Mandelstam variables are defined as 



•s 
t 
u 



(P1+P2) 2 
(pi - P3) 2 

(P2 - Ps) 2 ■ 

Following the simplified ACOT prescription ji?], HE we take m& = everywhere 
except in the Yukawa couplings. Then the tree level amplitudes of the s, t and u channels 
are: 

3C bh C hhh m b (fj,)Ml 



M = M?5 



s u Jt 



[s - Ml + iM h T h 



■Sjiv{p 2 )u{pi) 



M, 



v 2 t 
v 2 u 



where i and j are color indices for initial b and b quarks. The /i parameter is an arbitrary 
mass, which is introduced such that both the renormalized strong coupling and Yukawa 
coupling are dimensionless in N dimensions. The corresponding spin- and color-averaged 
matrix elements squared including interference terms are 

2 
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where Y h is the decay width of the Higgs boson. 

Summing the above terms, we obtain the total matrix element squared 



(|M | 2 ) 



<|M S °| 2 > + (|M t °| 2 > + (|M°| 2 > + 2{Re(M?M° u )) 
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r 2 
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The parton level cross section for inclusive bb 

1 1 



hh production becomes 
^LO = / ~(\M \ 2 )d$ 2 (bb ^ hh) 



(2) 



where d& 2 (bb — > hh) denotes the two-body phase space factor and a factor of 1/2 comes 
from identical particles in the final state. 



4 



III. NEXT-TO-LEADING ORDER CORRECTIONS FOR bb hh 

To determine the next-to-leading order (NLO) corrections for Higgs boson pair production 
in bottom quark fusion, we evaluate the cross section for (a) the leading-order subprocess 
(bb — > hh), (b) the a s corrections, and (c) the 1/A corrections (bg — > bhh), where A = 
\n(Mh/rrib). The order a s corrections have contributions from one-loop diagrams with virtual 
gluons (bb — ► hh) and tree-level real gluon emission (bb — > hhg). The NLO correction 
contains both the a s and the 1/A corrections. 

We write the parton level NLO cross section as 

^lo(xi,x 2 ,h) = a LO (xi,x 2 ,fi) +5a NL o(xi,x 2 ,fj) 
Sa NLO (x 1 ,x 2 , fi) = 5a as (x 1 , x 2 , fi) + 5a 1/A (x 1 ,x 2 , fx) (3) 

where <7lo(xx,x 2 , /i) is the leading order (Born) cross section and Sa^o(xi,x 2 , fi) is the 
next-to-leading order correction to the Born cross section, x 12 are the momentum fractions 
of the partons, X\X 2 = s/S, S is the center of mass energy of the hadrons and // = /zr is the 
renormalization scale. 

The 0(a s ) correction includes contributions from both virtual and real gluon emission: 

5a as (xi,x 2 , jj) = 8a v (x 1 ,x 2 ,fi) + 5a r (x 1 , x 2 , fi) (4) 

where 5& v (xi, x 2 , /i) and 8a r (xi, x 2 , //) represent virtual and real gluon emission NLO cor- 
rections to bb — > hh. 

A. Corrections with virtual gluons 

The one-loop diagrams for the 0(a s ) corrections to bb — > hh are shown in Fig. EJ The a s 
corrections involving virtual gluons are evaluated as the interferences between Born diagrams 
(Fig. and one-loop virtual diagrams (Fig. |2J). 

\M V \ 2 = M M loop + M loop M = 2Re(M loop M ) . 

We evaluate the amplitudes of the one-loop diagrams by applying dimensional regular- 
ization in N dimensions with N = 4 — 2e. The virtual diagrams are computed analytically 
and all tensor integrals are reduced to linear combinations of one-loop scalar functions. All 
amplitudes for the virtual corrections can be reduced to Born amplitudes multiplied by 
coefficients Xi. We find, 

2Re(M loop W ) = 2C F g 2 s Re{ 

where Cf = 4/3 and 

X s = Xg 

Xt = X\ + X3 + X5 + Xj 

X u = X 2 + X4 + Xq + Xg 



X S \M»\ Z + X t \M?\ 2 + X U \M" U \ Z + (X t + X u )Re(M?M°) } (5) 
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FIG. 2: One-loop virtual corrections to bb — > hh. 



The coefficients are, 
X 1 = -A (1 + 1- ln(-O) 
X 2 = -A (1 + 1- ln(-u)) 
X 3 = X 5 = 2A|(-t)- £ Q + 2 



X, X 6 2A{( //)'(;• • 



M, 2 



M, 2 



1 1 

X 7 = 2A -- + - 



X 8 = 2a{-4 + - 



2M 2 / -f \ . ' 

772 +ln(«) 



M 2 h -t \M 2 h/ 

2M 2 , /-u\ , . . 
— 5— — In — 7 + n(s 

^-;ln( S ) + iln 2 ( S ) + l-^ 



Mi, 



lln 2 (M 2 )-lln 



where A is a normalization factor 



r(l + e)(47r/i 2 ) e . 



16tt 2 
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The function F(x) is denned as 
—x 



F(x) 



Ml-x 



X 

sfih 




- 2 In 2 


f-x^/s 


\ Ml 



-ln 2 (M 2 ) + ln 2 (-s) +7T 2 
1-Ph 



In 



1 + Ph 



+ 2Li 2 



Ph 



2LU 



x 



2tt 2 



(7) 



where Ph = y 1 — AMl/s and is the dilogarithm or the Spence function 

The virtual corrections contain both ultraviolet (UV) and infrared (IR) divergences. In 



the MS scheme, the b quark Yukawa coupling is renormalized with the counter term |31l 



5mh 



-A- 



167TOJ, 



m b e 

This counter term contributes to the total matrix element squared by 



IM 



5m 



CT 



6 |M°| 2 + 4: 



32 Ana, 



, Srrib 
m b 



M t °| 2 + |M^| 2 + 2Re{M?M°) 



r 0|2 



M°| 2 + 2(|M t u | 2 + |M^| 2 + 2ReM?M%) 



r0|2 



1 2 



Summing over all relevant contributions, we obtain the following expression for the one- 
loop virtual corrections 



IM, 



2Re(M loop M ) + |M CT | 



M \ 2 U 



647TOJ. 



where IM^I 2 contains the finite terms: 



|M, 



D\ 



\n 2 (s) + 



2tt 2 



- 1 



im: 



1 1 W A 3 

? + 7 ln(s) "27 



1 2 



+ A- 



QAira. 



IM; 



D\ 



(8) 
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+ 



Ml 



Ml-t 

Ml 
Ml-u 



\n\Ml) + ln 2 (-t) +7T 2 |M°| 2 + 



ln 2 (M 2 ) + ln 2 (- M ) + 7r 5 



IMT + 



Fit) + \-\ H-t) 

F(u) + ^-^ln(-u) 



IM 



0|2 



IM' 



1 2 



Mi 



+ 



(Ml-t)(Ml-u) 
-ln(-t) - ^ ln(-w) + 2tt 2 + 3 
Ml 



ln 2 (-t) + ln 2 (-u) + F(t) + 



(M 2_ t)(M 2^){-™-ln 2 H)-tln 2 (-.) 

(t + «)Ue(M t M0) 



F(t) + F(u) - ^ ln(-t) - ^ ln(-tt) + tt 2 + 3 



+ 



tu 



F(t)+F(«) + 3-|lQH)-|lB(-u) 



Re(M?M°) .(9) 



(M?-t)(M|- W ) 

The divergences left in Eq. |H] are infrared divergences which are canceled by the infrared 
divergences in the real gluon emission corrections discussed in the next subsection. 
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B. Real gluon emission 



The Feynman diagrams for real gluon emission (bb —>■ hhg) are shown in Fig. |3J We 
assign the momentum as 

KPi)KP2) -> h(p 3 )h(p 4 )g(p g ). 

The real gluon emission corrections have infrared (IR) and collinear singularities. The 
infrared singularities cancel the virtual infrared singularities in Eq. |H] and the collinear sin- 
gularities are absorbed into parton distribution functions. We apply the two cut-off phase 
space slicing (PSS) method |2l|, |22| to isolate these singularities in different regions of phase 
space. 

We introduce two small cutoffs to split the phase space in the real gluon emission process. 
First, we introduce a soft cutoff (S s ) to separate the phase space of the process bb — > hhg 
into soft and hard regions according to the emitted gluon energy. The soft region is the 
region where the radiated gluon energy satisfies E g < <5 S ^, while the hard region is the 

region where the gluon energy satisfies, E g > 8 S ^ . Then the contributions from real gluon 
emission can be decomposed as 

&r = Psoft + Phard ■ (10) 

Although the energy of the emitted gluon in the hard region is above the threshold, there 
still exist singularities when the emitted gluon is parallel to one of the initial bottom quarks. 
We introduce a collinear cutoff (5 C ) to isolate this collinear singularity. The phase space in 
the hard region is decomposed into hard/collinear and hard/non-collinear regions. In the 
hard/collinear portion, the gluon is emitted within an angle satisfying 

— — < d c or — — < b c . (11) 

The parton level cross section in the hard region is split into hard/collinear and 
hard/noncollinear regions, 

&hard = &hard/c ~\~ &hard/nc j (12) 

where ahard/c is obtained by integrating over the hard/collinear region of the gluon phase 
space. It includes the collinear singularities and can be evaluated analytically in N dimen- 
sions. The hard non-collinear cross section (d-h ar d/nc) is finite. We have calculated crh ar d/c 
and ahard/nc numerically with a collinear cutoff (5 C ). The analytic matrix elements squared 
are used in our calculations for a s and 1/A corrections. Since 8 S and S c are arbitrary cut-offs, 
the dependence of the cross section on S s and S c is not physical and cancels in the total NLO 
cross section. 



1. Soft gluons 

The soft gluon emission amplitudes for the diagrams in Fig. |H1 are obtained by setting 
the gluon momentum p g to zero everywhere except in the denominators that are singular as 
p g —>■ 0. The soft gluon emission corrections are 

Mso f t = 9 2 X 3l - (M s ° + M,° + Ml) . (13) 
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FIG. 3: Feynman diagrams for Higgs pair production in bottom quark fusion with real gluon 
emission (bb — > hhg). 



There is a logarithmic divergence in the integral of the matrix element squared over 
the soft three-body phase space. The three-body phase space in the soft gluon emission 
approximation is, 



3 I soft 



2E 3 (2tt) n ' 1 2E 4 (2tt) n - 

(d$ 2 )(d$ fl U /t ) 



l {2Ti) N 5 N (p 1 +p 2 -p,- Pi ) 



(14) 
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where d$2 is the two-body phase space factor and d§ g \ so ft is the soft gluon phase space. In 
the center of mass frame of the incoming partons, 



9 1 soft 



r(l-e) 7r e 
r(l -2e) (2tt) j 



2 ° s 



dE g E]-* 



d6 1 sin L - 2€ 6 1 / dQ 2 sirr le Q 2 . (15) 



Together with the matrix element in the corresponding soft approximation, this integral can 
be evaluated analytically, 



IM: 



/ |2 



) 



g,| a0 /i(|M so /t 



2 ) 



4na s (\M \ 2 )C F 



1 /4vr/i 2 V r(l-e) 



|M | 2 )A 



647ra, 



4vr 2 e 2 V / T(l - 2e) 
3--^(0--^) + Jln 2 K 2 )-^ 



(16) 



The divergences in Eq. [H)] cancel the IR singularities in Eq. |S] of the virtual corrections. 
Adding Eq. |S] and Eq. UHl together we obtain 



(|M„| 2 ) + (|ilC /t | 2 > = (|M | ) < A 



6 Air a ^ 



H® + 



+A- 



7i- 

y 



+ A- 



IM, 



(17) 



However, this equation still has a collinear singularity, which can be absorbed into the parton 
distribution functions. We discuss this collinear singularity in the next subsection. 



2. Hard gluons 



In the hard/collinear region, the hard gluon is emitted collinearly to one of the initial 
partons. The phase space is greatly simplified in the collinear limit. The initial-state b quark 
splits into a hard parton b' and a collinear hard gluon g by b — > b'g with approximately 
p b > = zp b and p g — (1 — z)pb. In the hard/collinear limit, the matrix element squared for 
bb — > hhg factorizes into the Born matrix element squared and the Altarelli-Parisi splitting 
function for b — > b'g 



(E \M hard /c\ 2 )(bb -> hhg) -> (4vra s )(£ |M | 2 ) ^ (/i 2 )^ + (1 <- 2) 

where is the Altarelli-Parisi splitting function for b — > at the lowest order 

'1 + z 2 



(18) 



-FfM-z, e) 



el 



1-z 

P bb (z) + ePl b (z) 



(19) 
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The process can be factorized in two steps. First, one incoming b quark radiates a hard 
gluon and becomes b'. Then this b' collides with another incoming b quark to produce two 
Higgs bosons by b'b — > hh. The hard/collinear phase space is, 



d<t> 



3 \hard/c~ 



d^ 2 (zp! + p 2 -> p 3 + Pa) - 



d^p 9 



(20) 



To carry out the integration over d N ~ 1 p g in the collinear approximation, we introduce a new 
variable, s bg = 2pi ■ p g , which is constrained by < s bg < |(1 — z)S c . The hard/collinear 
phase space of the gluon becomes, 



d 



N-l„ 



I hard I c 



(47T) 



1 



16vr 2 r(l - e 



-dzds bg [(l - z)s bg 



(21) 



The dsbg integral can be evaluated to find the cross section in the hard/collinear region, 



&hard/c 



1 



X 



dx\dx 2 b(x2 
1-z) 2 



as r(i- 



27rr(l -2e)\ s 
+ (b «-> 6) + (1 <-> 2) . 



e 



— b(x l /z)P bb (z,e)a h0 

xi Z 



(22) 



This equation has collinear divergences which we remove by absorbing them into bare 
parton distribution functions. We introduce a modified (NLO) parton distribution function 
at the factorization scale fip in the MS scheme: 



2n 

b(r) = ?H.r. / , / .)<|l + ^(47r)T(l + e) 



Ml 



27rr(l-2er ; Ve 



1 58 P bb (z)-b(x/z) . 



(23) 



Replacing b(x) in the lowest order hadronic cross section by b(x,fip) and dropping terms 
higher than 0(a s ), we obtain the Born cross section that is proportional to a s 



0"Bc 



dx 1 dx 2 b(x 1 , fi F )b(x 2 , ^f)^lo(x 1 ,x 2 , fi R ) 



+ / dxxdx^X!, fi F )b(x 2 , [i F )a LO (x 1 ,x 2 , fi R ) 



4« f 
3tt 



+ 



x 



dx x dx 2 b{x 2 , n F )a LO (x 1 ,x 2 , fx R 



a, ra 



27rr(l - 2e) 



(4tt)T(1 + e) 



(Any - 



H* 2 s) + 



i— $ s dz - 
Pbb(z) — b(x l /z,/i F ) + (b^b) 



+ (1 «-> 2) 



(24) 



The 1/e poles in the second line cancel the collinear singularities of the soft gluon corrections 
in Eq. El while the divergences in the third line cancel with the hard/collinear divergences 
in Eq. 1221 These cancellations leave a finite cross section that has dependence on /j, R and 
fi F . 

The remaining region has hard/non-collinear gluons and yields a finite contribution to 
the cross section. We compute this contribution numerically. 
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C. Corrections from bg — ► bhh 



Now let us consider the contributions from the parton subprocess bg — > bhh, which is an 
0(1/ A) correction to the LO process of bb — ► hh. In calculating this cross section, again 
we ignore the bottom quark mass, raj, except in the Yukawa couplings where the running 
mass is used. There are no IR singularities in the bg — ► 6/t/i subprocess. However, when 
we integrate over the momenta of the b quarks, there are initial state collinear singularities 
which arise from gluon splitting to a pair of collinear b quarks. These singularities are 
absorbed into gluon parton distribution functions. 




(1) (2) 




(5) (6) 




FIG. 4: The lowest order Feynman diagrams for bg — * bhh. 



The eight diagrams for the bg — >• bhh subprocess are shown in Fig. |3J We assign momenta 
to partons as 

b(pi)g(p2) -> b(p b )h(p 3 )h(p 4 ) . 
12 



As there are no IR divergences, we do not need to separate the final b quark phase space 
into soft and hard regions. The cross section is: 



(Jbg — J dxidx 2 b{xi)g{x 2 )a{bg — > bhh) + (1^2) 



(25) 



As shown in diagrams (1), (2) and (7) of Fig.^J the initial gluon can be considered to split 
into two b quarks. When these two b quarks are parallel to each other, collinear singularities 
appear. To remove these collinear singularities, we introduce a collinear cutoff (S c ) to split 
the final b quark phase space into collinear and non-collinear regions. In the collinear region, 
the final b is emitted within an angle satisfying, 



~(P2 -Pbf 



< 5 r 



(26) 



Using the same method as for the bb — > hhg real gluon emission corrections, we find the 
cross section in the collinear region: 



2vrr(l - 2e) 



x / — g(x 2 /z)P bg (z,e) 

1X2 Z 



2z 



(27) 



Pb g is the Altarelli-Parisi splitting function for g — > bb 1 at lowest order, 



p b 9 (z,e) 



z 2 + (l- z? 



1 r 

2 

P b9 (z) + ePl g (z) 



ez(l-z) 



Again we introduce a modified parton distribution function, 

/ \ ot s r(l — e) . / 1\ r 1 , .dz , , , 
g(x,fi F ) = ff(*) + ^ r(1 _ 2e) W {—J I P bg (z)-g(x/z) 



(28) 



(29) 



The contribution to the NLO total cross section from the bg initial state is then, 

en f dz 

- 1 / dxidx 2 / — crLob{xi,fi F )g(x 2 /z,fi F ) 

Z7T J Jx 2 -2 



+ (1 - -2) 



hi 



^ z 2 



+ J dx\dx 2 b(xi, fj, F )g(x2, fj, F )a ri 



z(l-z) 
> 6/1/1) + 



1 



2) 



(30) 



Here &lq is the Born cross section for bb — > and <j„ c is the Born cross section for 5g — > 
The 6g — > process has exactly the same contribution as bg — > 
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D. Next-to-leading order cross section for bb — ► hh 



We have obtained the virtual and the real gluon emission 0(a s ) corrections to bb — > hh. 
The real gluon emission corrections include three regions: soft, hard/collinear and hard/non- 
collinear. The total next-to-leading order cross section can now be assembled from the lowest 
order cross section, the virtual corrections, and the real gluon emission corrections. Summing 
all of the above pieces, we obtain the 0(a s ) next-to-leading-order cross section corrections, 



C"Born + <5cr Qs — 0Born + &v + &soft + ^hard/c + &hard/r, 



dx\dx2 b(xi, fx)b(x2, //) \ &lo 



4a« 
3tt 



ln(/i 2 ) (ln(5 s 2 ) + ~) 



& finite 



a s f ri—fis dz r - - 

+—C F dx x dx 2 — a LO HxJz, fi)b(x 2 , fi) + b(x 2 , fi)b(x 1 /z, fi) 
2tt J Jxi z L 



X 



1 + z 2 



111 



+ 1 



l-z \fi 2 z 2 
+ J dx x dx 2 b(x 1 , fi)b(x 2 , n)a hard/nc (bb -> hhg) 
+ (1 - 2) 



(31) 



where o finite represents the finite cross section with contributions from both virtual and soft 
gluon radiative corrections, 



"finite 



1 14a, 



2s 2 3tt 



1 7T 2 

2 v s! 3 



|M | 2 + \M D \ 2 \d<& 2 {bb -»• /i/i) 



(32) 



We have taken /i = /x^ = //p. 

The C?(l/A) corrections from the subprocesses 6^ — > and 6(7 — > include collinear 
and non-collinear contributions and yield the 1/A corrections for Higgs pair production 
associated with one b quark. 



5a 1/A = 5a b ^ LO + 6a b ^ LO 



(33) 



Summing Eq. |^ and 1221 we get the total next-to- leading-order cross section for Higgs pair 
production from bottom quark fusion, 



Onlo = (y Born + 5a aa + Sax 



/A 



(34) 



The associated Higgs boson pair production with bb occurs via gluon fusion gg — ► bbhh 
and quark- ant iquark annihilation qq — > bbhh. These are subleading corrections to the NLO 
results given above. To estimate the cross section from these subprocesses, we have applied 
a nonzero value of m;, = 4.7 GeV for the bottom quark mass except in the Yukawa couplings 
where the running mass is used. We use MadGraph [2^ and HELAS |24| to calculate the 
cross sections for these processes and find that the contribution from qq — > bbhh is much 
smaller than that from gg — > bbhh. 
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IV. RESULTS FOR HIGGS PAIR PRODUCTION IN BOTTOM QUARK FUSION 



In this section, we present our results for the next-to-leading-order inclusive cross section 
for pp —>■ hh + X via bottom quark fusion, bb — > hh. We use the lowest order CTEQ6L1 
parton distribution functions (PDFs) |32j at the factorization scale fip with the leading- 
order evolution of the strong coupling o: s (//r) at the renormalization scale [ir to calculate 
the LO cross section and the CTEQ6M PDFs at /ip with the next-to-leading-order evolution 
of a s (/i#) to evaluate the NLO inclusive cross section. 

6 C = 5 s /^0, jli r = /li f = M h /2 
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FIG. 5: Order a s corrections Sa as (bb — > hh) (dash-dot, blue) versus the soft cutoff 5 S with \/S = 
14 TeV, hr = n F = M h /2 and S c = 5 S /W for (a) M h = 120 GeV and (b) M h = 200 GeV. These 
corrections have contributions from virtual gluons (cr v ), soft gluon emission (cr s ), hard/collinear 
gluon emission (cr/j arrf/ / c ), and hard/non-collinear gluon emission (cr hard / nc ). These graphs show the 
cancellation of the S s dependence between a v + a s + cr hard / c (dash, magenta) and cr hard / nc (dot, 
green). Also shown is the sum of <TB or n and 5a as (solid, red). 



We have introduced two arbitrary small cutoffs, the soft cutoff S s and the collinear cutoff 
5 C , to split the phase space in the real gluon emission corrections into soft, hard/collinear 
and hard/non-collinear regions. These separations are not physical and our final results are 
reliable only if they are not sensitive to these parameters. 

To check the dependence of the 0(a s ) NLO cross section on these two parameters, we 
present the order a s corrections for Higgs pair production 5a aa (bb — > hh) versus 6 S in Fig. 
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FIG. 6: Order a s corrections 5a a 
8 C with y/S = 14 TeV, [ir = = 



(bb — * hh) (dash-dot, blue) versus the hard/collinear cutoff 
M h /2 and 5 S = W5 C for (a) M h = 120 GeV and (b) M h = 



200 GeV. These corrections have contributions from virtual gluons (cr v ), soft gluon emission (<r s ), 
hard/collinear gluon emission (cr hard / c ), and hard/non-collinear gluon emission (cr hard / nc ). These 
graphs show the cancellation of the S c dependence between a v + a s + cr hard / c (dash, magenta) and 
a hard/nc (dot, green). Also shown is the sum of OBorn and 5a as (solid, red). 

and Sa as (bb -> hh) versus S c in Fig. El with = 14 TeV for (a) M h = 120 GeV and (b) 
Mh = 200 GeV. These corrections have contributions from virtual gluons (cr v ), soft gluon 
emission (<7 S ), hard/collinear gluon emission ((Jhard/c), and hard/non-collinear gluon emission 
(chard/nc)- These graphs show the cancellation of the 5 S dependence and the 5 C dependence 
between a v + a s + <7hard/c and crhard/nc- Also shown is the sum of <TBorn and 5a as . We have 
chosen S c = 5 S /10 and the renormalization/factorization scales to be /j,p = = Mh/2. 
There are several points to note. 

• We divide the phase space of the real gluon emission correction into soft, hard/collinear 
and hard/non-collinear regions by introducing two small cut-offs 5 S and 8 C . The cor- 
rection in each region is then integrated over the corresponding phase space. As we 
can see in these graphs, the corrections in each region are very sensitive to the values 
of S s and 5 C . 



Since these two small cutoffs, 5 S and S c , are arbitrary, the total cross section should 
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not depend on either one of them. These two figures show the cancellation of the 5 S 
and S c dependences. 

• At the renormalization/f actor ization scale fiR = fiF = Mh/2, the 0(a s ) NLO cross 
section corrections are comparable with the Born cross section. This implies that the 
0(a s ) corrections significantly increase the LO cross section. 
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FIG. 7: Order 1/A corrections 5a(l/A)(bg —* bhh) (dash-dot, blue) versus the collinear cutoff 
5 C with v 7 ^ = 14 TeV and hr = (m f = M h /2 for (a) M h = 120 GeV and (b) M h = 200 GeV. 
This figure shows the cancellation of the 5 C dependence between a c (dash, magenta) and a nc (dot, 
green). Also shown is the sum of OBorn and (solid, red). 

Figure E| shows the independence of order 1/A corrections (£<7i/a) on the collinear cutoff 
S c in the bg — > bhh subprocess. Like the a s corrections to bb — > hh, the contributions from 
collinear and non-collinear regions are very sensitive to the values of 5 C , but the total 0(1/ A) 
NLO cross section correction is independent of 5 C . The 0(1/ A) contributions from bg —* bhh 
are much smaller than the 0(a s ) NLO corrections from the bb initial state. 

In Fig. |S1 we study the dependence of the LO and NLO cross sections on the renor- 
malization and factorization scales. We present the next-todeading order cross section 
0"nlo(P£> — > hh + X) via bottom quark fusion versus /i = Hr = jiF with 5 S = 1CT 3 
and 5 C = 1(T 4 for (a) M h = 120 GeV and (b) M h = 200 GeV. Also shown are the LO 
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FIG. 8: Next-to-leading order cross section cfnlo(pP hh + X) (solid, black) in bottom quark 
fusion versus renormalization/factorization scale fi = fiR = hf with V~S = 14 TeV, 5 S = 10~ 3 and 
5 C = 1CT 4 , for (a) M h = 120 GeV and (b) M h = 200 GeV. Also shown are the LO cross section cr LO 
(dash-dot-dot, green), a s corrections Sa as (dash-dot, blue), 1/A corrections 5<r 1 / A (dot, magenta), 
and the NLO correction <5<7nlo = ^q s + <5<7i/a (dash, red). 



cross section (olo), a s corrections (5a aa 

<5o"NLO = 5<J as + &71/A- 

We note that: 



1/A corrections (Sai/\), and the NLO correction 



• The total NLO cross section corrections decrease with fi. The larger the Higgs mass, 
the slower the decrease of the NLO corrections. 

• The NLO cross section has less dependence on the renormalization and factorization 
scales than does the LO cross section. 

Fig. El shows the NLO total cross section pp — > hh + X in bottom quark fusion as a 
function of the Higgs mass (M h ) at the LHC with = 14 TeV, 5 S = 10" 3 and 5 C = 10" 4 . 
The renormalization and factorization scales are chosen to be (a) fi = Mh and (b) fi = 
The NLO cross section correction for [i = is bigger than for \i = In addition, we 

present the LO cross section (olo), ol s corrections (5cr as ), 1/A corrections (5<Ti/a), and the 
NLO correction <5<7nlo = 5o" Qs + 5cti/a- 
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6 S = 10" 3 , 6 C = 10" 4 



(a) fx n = fi F = M h (b) M R = /x F = M h /2 
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FIG. 9: Next-to-leading order cross section cfnlo(pP hh + X) (solid, black) in bottom quark 
fusion versus the Higgs mass (M/,) with y/S = 14 TeV, 5 S = 1CT 3 and 5 C = 10~ 4 , for (a) (j,r = 
[ip = Mh and (b) Hr = = Mh/2. Also shown are the LO cross section cjlo (dash-dot-dot, 
green), a s corrections 5a as (dash-dot, blue), 1/A corrections 5<tim (dot, magenta), and the NLO 
correction <5ctnlo = i5a as + 5<Ji/k (dash, red). 

In Fig. we show the LO and the NLO cross section including the 0(a s ) contributions 
from the bb initial state and the (9(1/A) contribution from the bg initial state. We sepa- 
rately show the 0(1/ A 2 ) contribution via gg — > bbhh as a function of the Higgs mass (Mh) 
at the LHC with \fs = 14 TeV. As in Fig. El we present our results at two renormaliza- 
tion/factorization scales, (a) /i = M h and (b) fi = M h /2 with 5 S = 1CT 3 and 5 C = 10~ 4 . 

In calculating the sub-leading contribution from the subprocess gg — > bbhh, we have 
evaluated the cross section numerically for pp — > bbhh + X via gg — > bbhh with a finite 
quark mass to regulate the collinear singularity. We take the b quark mass to be m?, = 4.7 
GeV. 

As is shown in Fig. both the LO and NLO total cross section for the bb — > hh process 
at a renormalization/factorization scale Mh are bigger than the corresponding cross sections 
at the scale M h /2. 

In the Standard Model, gluon fusion is the dominant source for producing a pair of Higgs 
bosons through triangle and box diagrams containing top and bottom quark loops. It has 
been demonstrated that it might be possible to study the trilinear Higgs coupling at the 
LHC through the gluon fusion production mechanism of Higgs boson pairs jE IE 0, IE ■ I n 
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FIG. 10: Leading order cross section cjlo (dash, blue) and next-to-leading order cross section ctnlo 
(solid, red) for pp — > hh + X via bottom quark fusion as a function of Mh at the LHC with \/~S = 14 
TeV for (a) fi = Mh and (b) fi = Also shown is the C(l/A 2 ) contribution from gg — > bbhh. 



the minimal supersymmetric standard model, bottom quark fusion could offer great promise 
to study the trilinear couplings of Higgs bosons for large values of tan f3 through bb — > hh. 

In Figure ^2 we show the cross section for pp — > hh + X via gluon fusion (gg — ► hh) 
through top and bottom loop diagrams as a function of at LHC with v^S" = 14 TeV for 
(a) fi = Mh and (b) fi = Mh/2. In the SM this production mechanism is significantly larger 
than the bb fusion mechanism. 



V. CONCLUSIONS 



In this paper we present complete 0(a s ) and 0(1/ A) QCD corrections to the production 
of a pair of Higgs bosons via bottom quark fusion at the LHC. We introduce two arbitrary 
small cutoffs, 8 S and S c , to compute 0(a s ) NLO soft, hard/collinear and hard/non-collinear 
gluon emission corrections to bb — > hh. The total 0(a s ) NLO cross section corrections are 
independent of the values of 5 S and S c . 

Our results are not sensitive to the difference between renormalization and factorization 
scales and we use the same renormalization and factorization scales. The LO cross section 
is very sensitive to the factorization scale. The cross section for the bb hh process with 
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FIG. 11: Cross section of pp — ► hh + X via gluon fusion versus the Higgs mass [Mh) at LHC with 
y/S = 14 TeV for the contributions from the triangle diagrams (dash, blue) and the box diagrams 
(dot, green) respectively. Also shown is the total cross section (solid). The renormalization and 
factorization scales are chosen to be (a) = fj,p = and (b) = fj,p = M^/2. 

large factorization scale is larger than the corresponding cross section with a low factorization 
scale. 

The rate for double Higgs production in the SM is very small, although the NLO cor- 
rections we have computed significantly increase this rate. However, the rate for Higgs pair 
production will be enhanced in models with large couplings of the Higgs bosons to b quarks. 
Our results are of interest in attempts to measure the trilinear Higgs coupling in such models. 
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APPENDIX A: HEAVY QUARK RUNNING MASS 



The two-loop running mass of a heavy quark in the MS scheme has the expression |17J,|2j 



where 



m(fi) = m(/i ) 
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l + a 1 ^M + (a? + a 2 )/2(^) : 
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where Nf is the number of quark flavors with mass less than > /io), and ( is the Riemann 
zeta function with ( 3 = 1.2020569. 



The MS mass at next-leading order is 



M Q = m(M Q ) 1 + C F — 



a. 



71 



(A3) 



Here Mq is the heavy quark pole mass. The evolution of the strong coupling can be 
found in Ref. I2i 
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